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ABSTRACT 

Physical processes regulating star formation in satellite galaxies represent an area of ongoing 
research, but the projected nature of observed coordinates makes separating different popula- 
tions of satellites (with different processes at work) difficult. The orbital history of a satellite 
galaxy leads to its present-day phase space coordinates; we can also work backwards and use 
these coordinates to statistically infer information about the orbital history. We use merger 
trees from the MultiDark Run 1 N-body simulation to compile a catalogue of the orbits of 
satellite haloes in cluster environments. We parametrize the orbital history by the time since 
crossing within 2.5 r v ; r of the cluster centre and use our catalogue to estimate the probability 
density over a range of this parameter given a set of present-day projected (i.e. observable) 
phase space coordinates. We show that different populations of satellite haloes, e.g. infalling, 
backsplash and virialized, occupy distinct regions of phase space, and semi-distinct regions of 
projected phase space. This will allow us to probabilistically determine the time since infall of 
a large sample of observed satellite galaxies, and ultimately to study the effect of orbital his- 
tory on star formation history (the topic of a future paper). We test the accuracy of our method 
and find that we can reliably recover this time within ±2.58 Gyr in 68 per cent of cases by 
using all available phase space coordinate information, compared to ±2.64 Gyr using only 
position coordinates and ±3.10 Gyr guessing 'blindly', i.e. using no coordinate information, 
but with knowledge of the overall distribution of infall times. In some regions of phase space, 
the accuracy of the infall time estimate improves to ±1.85 Gyr. Although we focus on time 
since infall, our method is easily generalizable to other orbital parameters (e.g. pericentric 
distance and time). 

Key words: galaxies: kinematics and dynamics, galaxies: clusters: general, galaxies: haloes 



1 INTRODUCTION 

We know that galaxies in clusters a re typically more 'red and dead' 
than t heir counterparts in the field dHogg et al. 2004; Balog h~et al.l 
2004), as wel l as dominantly elliptical (rather than spiral, see 
Dresslerlll980l) . This is thought to be due to a mechanism or com- 
bination of mechanisms that halt the collapse of cold gas into stars 
in a satellite galaxy as it orbits within a deep potential well - either 
by heating or removing the gas or by preventing the cooling of ad- 
ditional gas and consuming the existing supply. So me renowned 
mechanisms include ram pressure stripping (e.g. | Gunn & Gottl 
1 19721 : lAbadi et alJI 19991: Ijichvm et alj|2007t ISmith et alj|2010bl) . 
tidal stripping l Mayer et al.ll2006h . harassment dMoore et al Jl 19961 : 
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I Smith et al.l2010ah . st rangulation dLarson et al.ll9 80: Balog h et all 
hood) and mergers dToomre & Toomrel \\912i ICox et al.l 120061) . 
Internal quenc hing mechanisms have also been proposed - e.g . 
shock heating dBirnboim & Dekel|l2003l;lDekel & Birnbqiml 2006). 
AGN heating dCroton et al.ll2006l : lMcNamara et al.ll2006l) , see also 
Gabor et al. I d2010l) - which, while unable to account for the dif- 
ference between field and cluster galaxies, may each contribute to 
increasing the red fraction in both environments. 

It has so far been found difficult to produce a semi-analytic 
model for quenching that both reprodu ces the observed sta r forma- 
tion rate (SFR) d i stribu tion in clusters dWetzeletal.ld2012h . but see 
IWeinmann et alj j201Ch for some recent improvements). However, 
all the environmental quenching mechanisms listed above have at 
least one common characteristic : each is sensitive to the path taken 
through the cluster. Tidal stripping is strongest near the cluster cen- 
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tre. Ram pressure stripping varies with the local density of hot clus- 
ter gas and the relative velocity of the satellite, both of which vary 
radially in the cluster. Mergers are more probably in the outskirts 
of clusters. Harassment is most effective during high speed encoun- 
ters, which occur near the cluster core. Strangulation is triggered 
by the removal of the halo gas of the satellite, obsensibly via one 
of the aforementioned mechanisms. An environmental quenching 
model can therefore reasonably be expected to depend on the pre- 
vious positions and velocities, or more concisely the orbital history, 
of satellite galaxies. 

The orbital history of a galaxy is not directly observable; if we 
hope to compare environmental quenching models to observational 
data, we need some way to characterize the orbit of an observed 
galaxy. We might parametrize the orbital history of a galaxy in any 
of a number of ways - for instance by the time since infall into the 
cluster potential, the distance of closest approach to the cluster cen- 
tre or the number of orbits completed since infall - then attempt to 
find correlations between these parameters and observables using 
N-body simulation data. The observables we consider available are 
the distance between the satellite and the cluster centre projected 
on to the plane of the sky and the component of the velocity of the 
satellite relative to the cluster centre along the line of sight (LoS). In 
some nearby clusters it may also be possible to obtain distances to 
the cluster centre and the satellite via direct distance measurements 
(e.g. Tully-Fisher relation, surface brightness fluctuations, SNIa lu- 
minosities, etc.), but in general these will be much less precise than 
the other two coordinates - we will consider this information inac- 
cessible. 



Several author s dGao et al.l 12004 IWang et all 1201 it 
|Pe Lucia et alj |2012| ; iTaranu et al.l 1201 2|) have shown that the 
present radial distance from the cluster centre is negatively 
correlated t o the time s ince infall into the cluster and other orbital 
parameters. iGill et al l d2005h were amongst the first to compare 
the velocity distributi ons of different sat ellite populatio n s (see 
also IWang et "all l2005h , and more recently iMahaian et al.l d201 lh 
presented a method for separating satellite populations based on, 
amongst other parameters, their LoS velocities. 

This paper presents our method for reconstructing the 
parametrized orbital history of a satellite galaxy based on its 
present-day observable phase space coordinates; we defer its appli- 
cations to future papers. Here, we focus on one parameter - the time 
since infall on to a cluster-sized host - but stress that our method is 
easily adapted to other parameter choices. The time since the satel- 
lite first passed pericentre on its orbit in the cluster and the distance 
of closest approach to the cluster centre are two examples of other 
parameter choices that will be of interest in future applications of 
our method. 

In section [2] we describe the N-body simulation data and its 
conversion into halo catalogues and merger trees. In section [3] we 
present the method used convert the merger trees into our cata- 
logues of satellite orbits, and to a probability distribution function 
of times since infall given a pair of coordinates in phase space. 
In section|4]we discuss the importance of LoS velocity data in dis- 
criminating between different populations of satellites. We summa- 
rize in section [5] 

We assume the same cosmology used in the Bolshoi and 
Multidark Run 1 simulations with Q m = 0.27, Q \ = 0-73, 
fib = 0.0469, n s = 0.95, ho = 0.70, <r 8 = 0.82 jPrada et all 
12011 . 



2 SIMULATIONS 

To obtain a large sample of satellite orbits, we use the output of 
the MultiDark Run 1 (MDR1) dark matter only simulation with 
1 Gpch -1 box side length, 2048 3 particles, current cosmology 
(WMAP5/WMAP7), 8.7 x 10 9 /i _1 M mass resolution and 7 
kpc force resolution. The simulation was run from z — 65 to the 
present; the majority of snapshots are output at even intervals in 
scale factor a with some irregular intervals at small a. The resolu- 
tion in scale factor is of 0.0304 before a ~ 0.7 (z ~ 0.43) and 
doubles to 0.0152 afterwards (corresponding to a time resolution 
of about 0.210 Gyr at z — 0). This is taken into account in our 
discussion below. Full detail s of the simulation parameters are de- 
scribed in IPrada et al.l J2012I) . The sim ulation snapshots were pro- 
cessed by the ROCKSTAR halo fi nder dBehroozi et aL 201 lah and 
the merger tree code presented in lBehroozi et al.ld2011bl) . To pro- 
duce the data used in this paper, the merger tree code was slightly 
modified so that a halo may contain satellite haloes at distances of 
up to 2.5 times its virial radius (r\„ — r3eob) from its centre; the 
default code finds satellites within only 1.0 virial radii. This allows 
us to track sa t ellites out to th e largest apocen t ric distances which 
Mamon et al] d2004h (see also lGill et alj|2005l ; lBalogh et alj|200d ; 
Ludlow et al. 200§) show to be about ~ 2.5r2ooc(~ 2.2r360b) by 



using host-satellite linking as a proxy for cluster membership. 

In the following discussion, we will denote full 6D cluster- 
centric coordinates (r,v). The position and velocity centre of a halo 
is determined by averaging the positions and velocities of the sub- 
set of halo particles which minimizes the expected Poisson error in 
the coordinates, i.e. the particles occupying the region of highest lo- 
cal density (for more detail on how the coordinates are determined 
by the halo finder, we refer the interested reader to Behro ozi et al.l 
l2011al. §3.5.1). 

Projected coordinates will be denoted (R, V). Projection is 
done along the arbitrarily chosen third (z-)axis of the simulation 
box and includes a correction to the velocities to account for the 
Hubble flow. With this correction, the simulated velocity differ- 
ences can be directly compared to observation data, where ve- 
locity differences would presumably be measured using a red- 
shift difference. The projected distance between two points is 
R12 = \J (r2,x — ri >x ) 2 + (r 2 , y - ri tV ) 2 and the relative veloc- 
ity of point 2 with respect to point 1 is V12 = \(v2,z ~ v i,z) + 
H(r2,z — ri, x )\. Note that V > 0; this encodes our assumption 
that the distances to the two points are not known precisely, so only 
the magnitude of their relative velocity can be determined. 

For ease of comparison between satellites of different hosts, 
all spatial coordinates are normalized to the viria l radius r VII of the 
host h alo, which is defined using the formula of Brya n & Normanl 
( 1998): the radius enclosing a region with an overdensity of 360 
times the background density at z = 0. For readers more accus- 
tomed to normalization by r2ooc an approximate conversion at 
z = is r ™" c ~ 1.13. All velocity coordinates are normalized 
to the mis velocity dispersion a of the host halo, measured in 3D. 



3 METHOD 

First, the orbital history of a satellite and its host at z = are de- 
termined. Then the infall time of the satellite into that host halo is 
defined as the earliest time at which the satellite's progenitor iden- 
tifies the z = host's progenitor as its host. A host-satellite pair 
may be identified as soon as the satellite is within 2.5 times the 
virial radius of the host; though other criteria must still be met, 
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Figure 1. These two figures provide a measure of the importance of resolution effects in our analysis; see section [5] and i]4.4l for discussion. Left panel: 
The mass function of all satellite haloes in our initial sample (dotted black line) and the mass function of satellite haloes that have experienced at least one 
pericentric passage (solid black line), where in both cases the masses are measured at the time of infall on to their host. The coloured lines separate the satellites 
that have experienced a pericentric passage into bins of pericentric distance, as labelled. The total number of satellites contributing to each mass function are 
labelled N . We impose a mass cut as indicated by the thick vertical line (Men = 10 Mq), yielding a sample which is as complete as reasonably possible 
above Mcut. Right panel: Number density as a function of satellite radial position for various mass ratios fi, where the mass of the satellite is measured at the 
time of infall. The host halos have masses > 10 14 Mq . Satellites with low mass ratios are underabundant toward the centre of the cluster environment. 



as outlined in lBehroozi et al.l l l2011bh . in practice these are usually 
met immediately when the satellite crosses within 2.5r v i r . Only the 
orbital history with respect to the final (z — 0) host is considered, 
and in cases where a satellite has a hierarchy of hosts, intermediate 
hosts are ignored and the satellite identifies the largest as its host; 
in other words, the present work ignores "group pre-processing". 
Hosts are selected to be "cluster-sized", which we define as a halo 
mass > 10 14 Mq. With the MDR1 dataset, this yields a catalogue 
of ~ 570, 000 satellite orbits belonging to ~ 24, 500 different 
hosts. 

To ensure a sample complete in satellite mass, and that we 
are minimally sensitive to artificial disruption of satellite haloes, 
we impose a mass cut at M cul = 10 11 ' 9 Mq, where the mass is 
measured a t the time of infall. Us ing the stellar-to-halo mass ratios 
outlined in Behr oozi et al. d2012h . this corresponds to a cut in stel- 
lar mass at ~ 10 10 ' 3 Mq. The left panel of Fig.Q]shows the mass 
function for all satellites, and for satellites that have experienced at 
least one pericentric passage in bins of proximity of pericentric pas- 
sage. Our mass cut is safely above the mass resolution limit of the 
simulation. As in essentially all cosmological simulations, MDR1 
halos experience artifical d i sruption in high de nsity environments 
dKitzbichler &~W hite 20081: iKlvpin etaj]|l999h . This is the reason 
for the underabundance of lower mass satellites that have experi- 
enced a close approach to a cluster centre (blue curve in Fig. Q~| 
left panel). From these mass functions, we estimate that less than 
20% of halos with masses above our mass cut and pericentric dis- 
tances in the range 0.0 — 0.25r v i r have been artificially disrupted. 
Most of the artificially disrupted satellites are in the mass range 
10 11,9 - 10 131 M . We estimate that these missing satellites ac- 
count for less than 4% of the total halo population above our mass 
cut. The right panel of Fig.[T]shows the number density of satellite 
halos as a function of radial position for bins of satellite mass rela- 
tive to the host halo. Observations constrain the slope of this power 
law relation to be betw een —1.7 and —1.5, regardless of mass ra- 
tio dTinker et al. I l2012l) . The slopes shown in Fig Q] are somewhat 
steeper at about —1.9 for mass ratio /1 = 10" — 10°. The rea- 



son for this steeper slope is not precisely known, but the higher 
resolution Bolshoi simulation exhibits the same slope of —1.9 so 
we surmise that it is not due to a resolution effect. The key feature 
that we wish to highlight is that the slope (and shape) of this rela- 
tionship is mass ratio dependent in our dataset, with satellites that 
are smaller relative to their host being less abundant closer to the 
centre of the host. The underabundances of satellites highlighted 
by each panel of Fig.Q]are due to the same population of satellites; 
those which have orbited to within < 0.25r v k of their host and have 
a mass < 10 13 ' 1 Mq. See ^4.41 for a discussion of the impact of 
these missing satellites on our results. 

We impose one final cut, removing satellites that existed for 
less than 3 simulation snapshots before falling into a cluster. This 
prevents haloes near the mass resolution limit of the simulation 
from appearing suddenly inside a cluster and being assigned mean- 
ingless infall times. The remaining 242,790 satellites were binned 
in 100 projected position bins in 0.0 < R/r Y i r < 2.5 and 100 pro- 
jected velocity bins in 0.0 < V/a < 2.0. The set of infall times 
in each bin was used to create a probability distribution function 
(PDF) of infall times for each bin. 



4 RESULTS 

Before applying our PDFs to modelling environmental quenching 
(which will be the focus of a future paper), we should have an un- 
derstanding of a few systematic effects inherent in the method pre- 
sented in section [3] We will discuss the impact of projecting the 
data in both the radial and velocity coordinates in i]4.1| In i]4.2| we 
will present the PDFs and discuss some of their features. In ij4.3l we 
discuss the effects of both host and satellite mass on the distribution 
of satellites in phase space. Finally, in i]4.4| we discuss the impact 
of resolution effects on our results. 
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4.1 Projection effects 

While a simulation provides accurate values for all six phase space 
coordinates of an object, a typical astronomical observation can 
only measure three. The right ascension and declination give two 
spatial coordinates (the distance is unknown), while comparing the 
spectra of two objects can give the difference in velocity between 
them, but only for motion in a direction along the LoS. Since we 
ultimately wish to infer properties of observed objects from their 
coordinates, we must restrict our knowledge of simulated objects 
to these same coordinates. This can be achieved by ignoring one 
of the spatial coordinates of the simulation box - in our case, the 
third - and considering only the velocity coordinate corresponding 
to the ignored spatial coordinate. Additionally, we include a correc- 
tion to the projected velocity to account for the Hubble flow. This 
transformation applied to the spatial coordinates of a point x and 
its velocity u relative to a reference point y and its velocity w can 
be expressed: 

r = (x-y) = (an - yi,x 2 -y 2 ,x 3 - y 3 ) 
R = (xi — X2 — yi) 

V = (u — w) = (itl — Wl, U2 — W2, Us — IU3) 

=► V = (|(« 3 - w 3 ) + H(x 3 - 2/3)|) 

We predict the effect of projection on the radial coordinate by 
considering a random uniform distribution of points on a spher- 
ical shell. For this distribution the relationship between the ac- 
tual (r) and projected (R) radial coordinates is characterized by 

<# > = f ± vf-fg- ~ 0.79 ± 0.22 (la scatter). In our sam- 
ple of satellite galaxies, we find that {— ) agrees with this pre- 
diction to within 1%, both in the mean and in the scatter. The 
projected radial coordinate tracks the 3D radial coordinate more 
closely and more consistently at larger projected radii; at R ~ r v i r , 
(f ) ~ 0.83 ± 0.17 while at R ~ 0.1r vir , <f ) ~ 0.51 ± 0.31. 
The observed velocity coordinate does not have as straightforward 
a relationship with the actual quantity of interest - the radial com- 
ponent of the velocity difference between two points - but the in- 
formation about two components is lost instead of one, so we ex- 
pect a much larger typical difference between projected velocity 
and true velocity. We also lose the sign of the one remaining com- 
ponent of velocity since, without knowledge of the distances to the 
two points, we cannot know whether the distance between them is 
increasing or decreasing. 

The left panel of Fig.[2]shows the distribution of satellite infall 
timefl as a function of radial distance from the host centre mea- 
sured at z — before projection. Satellites that have recent infall 
times (near the top of the diagram) are necessarily concentrated 
near the edge of our definition of the cluster at 2.5 r v i r ; they have 
not had enough time to move anywhere else. A typical crossing 
time for our sample of clusters is about 6-8 Gyr, so most satellites 
that fell in 3-4 Gyr ago are near the centre of the host. Note that 
this seemingly long crossing time is due to our definition of the 
edge of the cluster at 2.5 r v „; the typical time to cross from r v i r to 
pericentre and back to r v ; r is about 2 Gyr. The satellite also spends a 
significant amount of time outside the virial radius after this initial 
crossing; it takes a further ~1 Gyr for a typical satellite to reach 
apocentre after making its first outbound crossing of r v h. 



1 The infall times produced by the method of section [3] occur at discrete 
times - those of the simulation snapshots. For the purposes of visualization 
only, some scatter was added to the times. 



There is a spread in the time taken to reach pericentre and the 
radii of the pericentres caused by the variety of possible orbits and 
details of the host potentials. The first apocentre after infall typi- 
cally occurs after about 6 Gyr; the population of objects between 
first pericentre and apocentre is termed 'backsplash' . Satellites with 
infall times earlier than ~ 7 Gyr have less distinct features in their 
distribution due to the increasing impact of variations in orbital his- 
tory, but the majority are confined within ~ 1 r v h-; we call this pop- 
ulation 'virialized'. We note an overall decreasing number of satel- 
lites with increasing time since infall. Some satellites with early 
infall times are disrupted by tidal interactions and do not appear in 
our merger trees. In other cases two infalling satellites may merge 
and appear in the trees as a single satellite. Finally, some halos may 
not host a galaxy. These three effects mean that observations of 
satellite galaxies around a cluster do not correlate perfectly with 
the distribution of dark matter halos around that cluster; this needs 
to be taken into account when applying our method to a sample of 
observed galaxies, but does not impact the method itself. There is 
one other contribution to the decrease in number of satellites with 
increasing time since infall. Because our mass limit of 10 11 ' 9 MqIs 
relatively large, at early times haloes above this limit were some- 
what rarer. 

The right panel of Fig. [2] shows the effect of projection in the 
radial coordinates on the same distribution as in the left panel. Fea- 
tures are shifted somewhat to lower radii (consistently with our 
expectation of (— ) = j) and broadened by the scatter about 
this mean deformation. All the populations and features discussed 
above are still identifiable. One notable omission from this diagram 
is any foreground/background objects, which are common in obser- 
vational samples, that could be confused with the satellite popula- 
tions. 

There is one other important effect to consider when inter- 
preting Fig. [2] (and others involving infall times). As a host halo 
accretes mass, its virial radius grows (slowly, except in the case of 
major mergers). Because of this, an orbiting satellite may appear to 
move further in coordinates scaled to r v „ than it otherwise would as 
the coordinates grow around it. This does not have a large impact 
on the positions, but contributes some of the scatter in the radial co- 
ordinate. One might consider choosing some radius that is constant 
with time to scale each halo, but other choices, such as the virial 
radius at the last snapshot, also introduce similar effects. 

Next, we consider the effect of projection on the velocity co- 
ordinates. The upper left panel of Fig. [3] shows the distribution of 
satellite haloes in phase space at z = 0. A typical halo would, given 
enough time, progress from large radii and low velocities down to 
low radii and high negative velocities, then switch to high positive 
velocity as it passes pericentre. From there it follows a series of pro- 
gressively shrinking concentric semicircles or chevrons, jumping 
from negative to positive velocity at each pericentric passage (for 
a more in depth theoretical background, see lBertschingerlll985l es- 
pecially fig. 6 therein). This normal progression is well represented 
by the distribution of our halo sample, however the individual or- 
bital 'shells' are not visible since they overlap, and we only expect 
1-2 shells given the orbital timescales and ages of these systems. 
Also, close encounters redistribute some haloes off this idealized 
track. 

Based on this expected movement through phase space in 
time, satellite haloes with different infall times should occupy dif- 
ferent regions of phase space. This is shown in Fig. [3] where the 
phase space distribution of haloes is plotted for a variety of bins 
in infall time. Fig. [4] shows the same binned distributions, but in 
projected coordinates, simultaneously showing the effects of both 
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radial distance r/r vir at z=0 projected radial distance R/r^ at z=0 

Figure 2. Left panel: Abundance of satellites in bins of cluster infall time, defined as first inward crossing of 2.5 r v ; r , and radial distance from cluster centre 
at t = 0. In the interval t = to t ~ 3.5 the 'infalling' population is visible, composed of satellites that have not yet experienced a pericentre passage. From 
t ~ 3.5 to t ~ 6 the 'backsplash' population is visible, composed of satellites that have passed pericentre once and are approaching apocentre (t ~ 6). The 
'virialized' population (t > 6) are on a second or subsequent orbit. Gray horizontal lines illustrate the times of simulation snapshots. Right panel: same as left 
panel, but using radial distance from cluster centre projected along one axis of the simulation box (simulating coordinates accessible in observations). This 
widens the distribution of radii of satellites at a given infall time, but the same populations as in the left panel are visible. 



radial and velocity projections, but emphasizing the latter. Much of 
the structure visible in Fig. [3]is lost in projection, but haloes with 
different infall times still occupy different regions of phase space. 

4.2 Infall time PDFs 

Since our ultimate goal is to model quenching, perhaps the most di- 
rectly useful question we can ask of our dataset is, given a position 
in projected phase space, what is the distribution of possible infall 
times (or another parameter of interest), and what is the likelihood 
of each? This question is easily answered in a statistical manner by 
sampling all satellites in a small region of phase space and binning 
them by infall time. The result is shown in Fig.|5]for a selection of 
points in phase space, varying both the radial coordinate (left-right 
across the panels) and the velocity coordinate (up-down). We focus 
here on the trends with the velocity coordinate. The first trend we 
note is that at a given radius, satellites with higher ul s typically 
have slightly more recent infall times than their low tt s counter- 
parts. This is because a satellite with high ul s is more likely to 
have a high total speed and can penetrate deeper into the host po- 
tential in a given amount of time than a satellite with low speed. 

In some cases, the velocity coordinate allows us to discrimi- 
nate between different satellite populations. Consider the rightmost 
panels (R — 1.5) of Fig. [5] One peak in the distribution of in- 
fall times is clearly visible at all values of uloS, at a fa 0.85; it is 
made up of infalling satellites. A second peak at a ~ 0.6 is made 
up of backsplash satellites, and is only present for low values of 
vlos- The backspash galaxies have lower kinetic energy relative to 
the host potential than recent infalls due to a combination of mass 
accretion by the host (increasing a, therefore causing an apparent 
slowing of haloes) and dynamical friction. 

Given a sufficiently large sample of observed galaxies, the 
PDFs we have created will allow us to statistically assign an infall 
time to each one using all the available dynamical information and 
a meaningful uncertainty in our assignment. We have developed a 
code and a compact data table to this end which we are prepared to 
share on request. 



To test the ability of our PDFs to correctly assign an infall time 
to a satellite halo, we use our PDFs to estimate the infall time of 
each satellite in our sample. This estimate is then compared to the 
actual infall time of the satellite, which we know from tracking its 
orbit. The distribution of At = tmfaSi (actual) — ti„f a u(estimated) is 
plotted in Fig.[6] A At of represents a correct guess, so a stronger 
peak about At = represents a higher rate of success. As a quanti- 
tative measure of the strength of the peak, we calculate the standard 
deviation of the distribution. By using our PDFs based on knowl- 
edge of the (R, V) coordinates of each satellite, we assign the in- 
fall time correctly to within ±2.58 Gyr in 68 per cent of cases. For 
comparison, we also plot the distribution of At obtained if only the 
position coordinate R of the haloes is known, accurate to within 
±2.64 Gyr in 68 per cent of cases, and if none of the coordinates 
are known; in this case, we simply draw values from the distribu- 
tion of infall times of our entire sample at random, and find we are 
accurate within ±3.10 Gyr in 68 per cent of cases. Finally, we plot 
a curve showing the distribution of At for a particular subsample 
of satellites chosen by eye which have > — | ± 2. This sub- 
sample is designed to retain primarily infalling satellites (see Fig. 
|4), and in this case the assignment of infall times is even more re- 
liable. This demonstrates that if a particular population of satellites 
is of interest, it is often possible to choose a region of phase space 
that maximises the likelihood of correctly assigning fjjfaii. For this 
example, we correctly assign iinfaii within ±2.48 Gyr in 68 per cent 
of cases. In practice, this makes it possible to increase the purity of 
an observational sample of galaxies at the expense of completeness 
of the sample by limiting the region of phase space from which the 
sample is drawn. 

Fig. [7] shows the accuracy of the assignment of tmSail at the 68 
per cent confidence level for bins in (R, V) space. The PDFs give 
the most reliable results in those regions with both good accuracy 
and a large population of satellites. In the outskirts of clusters, an 
accuracy At of about 3 Gyr is sufficient to robustly separate in- 
falling and backsplash satellites. 

We conclude that using all observable coordinates of a satel- 
lite (R, V) offers an improvement in assigning t-misM correctly over 
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Figure 3. The upper left panel shows the phase space distribution of satellite haloes (no projection). The infalling population of haloes is particularly distinct, 
forming the long dark bar with v < 0. The backsplash population forms the upper and lower edges of the rest of the main distribution, while the virialized 
population fills in the centre. The other panels correspond to bins by satellite infall time, as labelled (most bins are 1 Gyr). Each panel shows the z = 
distribution of satellites in phase space. Different satellite populations occupy distinct regions of phase space; compare for instance the upper right panel 
showing mostly infalling satellites, the centre right panel showing mostly backsplash satellites and the lower right panel showing mostly virialized satellites. 
Each bin is also labelled with the number of satellites N contained in the bin. 



using only the position coordinate R comparable to, but somewhat 
less than, the improvement seen when using the position coordinate 
R rather than no knowledge of satellite coordinates at all. Careful 
selection of the region of phase space to be sampled can, in many 
cases, further increase the reliability in estimating Wall- 



4.3 Mass trends 

While haloes of different masses are expected to be self-similar in 
most ways, we still find some trends with both host halo mass and 
satellite halo mass. We study these trends by separating the haloes 
of interest - either hosts or satellites - crudely into a 'high mass' 
and a 'low mass' bin, then compare the distribution of satellites in 
phase space and infall time for the two bins. 

Fig.[8]shows a comparison of satellites around high- and low- 
mass hosts in infall time vs. radial position space (left panel), and 
in phase space (right panel). Low mass hosts are defined to have 
10 14 M < M < 10 14 ' 5 M Q while high mass hosts are in the 
range 10 14 ' 5 M Q < M < 10 15 M Q . The population of satellites 
around each type of host is normalized, then the two are subtracted 



in each infall time - radial distance bin and compared to the fiducial 
abundance of satellites in that bin. Cells with a reddish colouration 
are preferentially occupied by satellites around high mass hosts, 
while cells with a bluish colouration are preferentially occupied by 
satellites around low mass hosts. 

In the infalling population of satellites (t < 4 Gyr in Fig. [8] 
left panel) the satellites of high mass hosts appear to fall into their 
host somewhat faster than the satellites of low mass hosts; at a given 
time since infall, a satellite of a low mass host has a larger typical 
radial distance than a satellite of a low mass host. We propose two 
possible explanations for this effect. First, satellites of high mass 
hosts tend to have more radial orbits than satellites of low mass 
hosts (Wetzel 201lh . causing them to move radially inward more 
rapidly on average. Alternately, because the hosts are continually 
accreting mass, their virial radius gradually increases. Hosts with 
higher masses have higher present-day accretion rates; this is vis- 
ible in Fig. [8] (right panel) where the region of phase space oc- 
cupied by infalling satellite s shows an excess of satellites around 
higher mass haloes (see also lWechsler et al.l2002l fig. 3 therein es- 
pecially). Their satellites therefore appear to move inward slightly 
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Figure 4. As in Fig. [5] but using radial distance from cluster centre projected along one axis of the simulation box and velocity projected along the same axis 
(again simulating coordinates of an observed system). Much of the structure visible in Fig.[3]is lost in projection, but the use of LoS velocity data allows better 
spearation of different satellite populations than using the radial coordinate alone. The dashed gray line marks = — ^ + 2; see §4.2. 



faster than those of low mass hosts in this coordinate system. These 
two explanations could be considered an argument for a different 
normalization of the radial coordinate, however other choices (e.g. 
the virial radius of the host at the infall time of each satellite) still 
introduce similar effects, and the present virial radius of the host 
has the advantage of being more related to observable quantities. A 
more practical approach is to compute PDFs for narrow host mass 
bins (~ 0.5 dex) and use whichever is applicable to an observed 
system of interest; our sample contains enough satellites to make 
this feasible. 

We also consider trends with satellite mass. Because satellite 
mass typically decreases with increased time spent in a cluster envi- 
ronment, we use the satellite mass at time of infall (first crossing of 
2.5 r v ir) as a characteristic mass for this analysis. Fig. |9]shows that 
lower mass satellites experience backsplashes to larger radii than 
high mass satellites. This is because the slowing due to dynamical 
friction is proportional to the mass of the satellite, so that higher 
mass satellites lose a larger fraction of their kinetic energy during 
a passage through a cluster. To account for this trend with satellite 
mass, separate PDFs can be produced for relatively narrow bins in 
satellite mass. 

Care must be taken when using our method to interpret the 



distribution of galaxies in a cluster. Baryons (especially the stellar 
component) are typically more tightly bound than their associated 
dark matter halo, and should therefore outlive the halo in the cluster 
environment. This restricts the satellites, as traced by galaxies, that 
can be studied using our method to those which we expect to have 
a surviving DM halo. With the MDR1 dataset, this corresponds 
to galaxies with associated DM halo masses of~ 10 119 M Q or 
greater. Smaller satellites may still be studied by usin g higher res- 
olution simulations (e.g. the Bolshoi simulation of iKlypin et al.l 
(2011), which has identical parameters to MDR1 but with higher 
resolution and a smaller box size). Low mass satellites are much 
more abundant, so the smaller volume of a higher resolution simu- 
lation should not impact our ability to obtain a large sample, except 
perhaps for low mass satellites around very large hosts (of which 
there may only be a few in a simulation such as Bolshoi). 



4.4 Resolution Effects 

As shown in section [3] some halos which experience a close ap- 
proach to the centre of a larger halo experience artificial disruption. 
We estimate that any given bin in our PDFs is less than 20% in- 
complete, and that our sample of satellites as a whole is less than 
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Figure 5. Probability density functions (PDFs) of infall times for a selection of points (R/r vlI , V/tr) in projected phase space. Each point is sampled in a 
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Figure 8. Relative abundance of satellites in the space of infall times vs. radial positions (left panel) and phase space (right panel), comparing satellites around 
low mass (10 14 — 10 14,5 Mq) hosts and high mass (10 14,5 — 10 15 Mq) hosts. The two populations are normalized, subtracted in each bin, then compared 
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Figure 9. Relative abundance of satellites in the space of infall times vs. ra- 
dial positions comparing low mass (10 11,9 — 10 12 ' 2 Mq) and high mass 
(10 12,2 — 10 12,8 /i — 1 Mq) satellites. The two populations are normalized, 
subtracted in each bin, then compared to the fiducial abundance of satel- 
lites in that bin (arbitrarily chosen to be the high mass population). Only 
bins containing at least 100 satellites are shown. High mass satellites are 
slowed more than low mass satellites by dynamical friction and therefore 
have lower typical backsplash distances than their low mass counterparts. 



4% incomplete down to 10 1 Mq. The effect on the predictive 
power of our method was investigated by generating an orbit cat- 
alogue with extra orbits added to compensate for artificially dis- 
rupted halos. The effect of artificial disruption on the distributions 
of At (such as those shown in Fig.[6} is very small, causing changes 
on the order of a few tenths of a percent. 

Our initial interest for applications of our method is to large 
satellites (halo masses > 10 12 ' 5 Mq) around hosts with halo 
masses > 10 14 Mq, which motivated our choice of the MDR1 
simulation as a starting point. Applying our method to smaller 
satellites of smaller hosts is simply a matter of choosing an appro- 
priate simulation; our method is easily applied to any simulation 
that can be processed by the ROCKSTAR code. The Bolshoi simu- 
lation is of particular interest since it is identical to MDR1 in all 
respects save resolution and could therefore be used in conjunction 
with MDR1 to extend the range of our method to lower mass ranges 
while maintaining good statistics at the high mass end. 



5 CONCLUSIONS 

The phase space distribution of infalling, backsplash and virialized 
satellite halos are different but not everywhere distinct. The LOS 
v elocity distributio ns we recover are in agreement with the results 
of lGill et al.U2005h . Like them, we find that different populations of 
satellite halos are better separated in phase space than in the radial 
coordinate alone, but that there is no immediately obvious cut that 
we can impose on the projected phase space coordinates of a satel- 
lite galaxy to separate different populations. iMahaianetallfeOllh 
identified some regions of projected phase space where parts of the 
infalling or backsplash populations could be picked out with lit- 
tle contamination by other populations. Building on this idea, we 
have examined the entire projected phase space and determined the 
confidence with which the time since infall can be assigned to a 
satellite occupying that region. This lets us easily identify regions 



© 0000 RAS, MNRAS 000, 000-000 



10 K. A. Oman, M. J. Hudson & P. S. Behroozi 



where the time since infall can be determined accurately enough to 
reliably separate satellite population s. Such regions that w e identify 
are a superset of those identified bv lMahajan et al 1 d201lh . This al- 
lows an increase of the area of projected phase space and therefore 
the number of satellite galaxies available for use in studying star 
formation and other effects that may depend on the orbital history 
of a satellite. 

We have developed a tool that estimates the infall time and 
a confidence in the estimate for a satellite halo given its projected 
phase space coordinates. The possible infall times given a pair of 
projected phase space coordinates are weighted according to the 
frequency of their occurence in simulation, then an infall time is 
randomly selected from this weighted distribution. We predict that, 
when applied to a large sample of observed galaxies, this method 
will allow correlations between infall times and satellite star forma- 
tion histories to be studied. Our method is easily adapted to other 
similar parameters (closest approach to host centre, time since first 
pericentre, etc.) that will be considered in our forthcoming study of 
SF quenching in cluster environments, and to 'preprocessing' sce- 
narios - accretion on to a group sized halo before accretion on to a 
cluster sized one. It would be straightforward to adapt our method 
to systems at higher redshift where we think the different satel- 
lite populations may be better separated than at low redshift. Our 
framework would also be very well suited to studying dark matter 
stripping of satellites. 
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